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TOh to* 


BASIC NOTATION 


m = the mass of the disturbed body, 

M = the mass of the Sun, 
f = the gravitational constant, 

= f ( M + m) , 

r = the heliocentric position vector of the disturbed body, 
r = |r|, 

r° = the unit vector along r, 

n° = the unit vector normal to r and lying in the orbital plane of the disturbed body, 
a = the semi-major axis of the orbit of the disturbed body, 
e = the eccentricity of the orbit of the disturbed body, 
g = the mean anomaly of the disturbed body, 
e = the eccentric anomaly of the disturbed body, 

p = a(l “c 2 ), 

j = the unit vector directed from the Sun toward the perihelion of the disturbed body, 
2 = the unit vector normal to P l and lying in the orbital plane of the disturbed body, 



A = the true orbital longitude of the disturbed body, reckoned from the departure point 
of the ideal system of coordinates, 

x = the true orbital longitude of the perihelion of the disturbed body in the ideal system 
of coordinates reckoned from the departure point, 


a = the angular distance of the ascending node from the departure point, 

Ri, r 2 , R 3 = the unit vectors along the axes of the ideal system of coordinates. r 2 and r 2 are in 
the osculating orbital plane of the disturbed body, R 3 is normal to this plane. The 
intersection of R : with the celestial sphere is the departure point. 


^3 ^1 X 


s,. s 2 ,s 3 


the initial values of R 1 ,R 2 ,R 3 , respectively. 

the Gibbs' vector. This vector defines the rotation of the orbital plane of the dis- 
turbed body from its initial position to the position at the given time t, 

the mass of the disturbing body, 

the heliocentric position vector of the disturbing body, 
the semi-major axis of the orbit of the disturbing body, 


v 


min am 


e' = the eccentricity of the orbit of the disturbing body, 
g' = the mean anomaly of the disturbing body, 
e # = the eccentric anomaly of the disturbing body, 

P x ' = the unit vector directed from the Sun toward the perihelion of the disturbing body, 
P 2 ' = the unit vector normal to p x ' and lying in the orbital plane of the disturbing body, 

A i = 
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A DISCUSSION OF HILL’S METHOD OF SECULAR PERTURBATIONS AND 
ITS APPLICATION TO THE DETERMINATION OF THE ZERO-RANK EFFECTS 
IN NON-SINGULAR VECTORIAL ELEMENTS OF A PLANETARY MOTION 


by 

Peter Musen 

Goddard Space Flight Center 


INTRODUCTION 

In the motion of planets and comets the purely M secular 1 ’ part of the disturbing function pro- 
duces perturbations in the elements with periods of many thousands of years. Poincare (1905) 
classified them as of ’’zero-rank.” They constitute a most essential part of those perturbations 
which regulate behavior of the orbit over an interval of, say, some millions of years. Knowledge 
of these perturbations of zero-rank is also valuable in cosmogony and paleoclimatology. Hira- 
yama’s discovery (1923) of families of minor planets represents one of the most beautiful cosmo- 
gonical applications of the theory of zero-rank secular perturbations in its linearized form. The 
paleo-climatological significance of the zero-rank effects in the motion of the Earth is explained 
in the work of Milankovitch (1948). 

In the case of artificial satellites the time scale becomes contracted, and the periods of the 
perturbations become only a few years instead of thousands. Shute (1964) noticed large oscilla- 
tions in the orbital eccentricity and inclination of satellites launched deep into cislunar space. 
Thus, knowledge of the zero-rank perturbations facilitates planning the launchings of satellites 
into elongated ellipses in cislunar space with lifetime prolonged or shortened as needed. 

In earlier times, beginning with Lagrange, information about the secular behavior of the orbital 
elements was obtained on the basis of linearized differential equations and under the assumption 
that the orbital eccentricities and the inclinations of both the disturbed and disturbing bodies re- 
main small. Since then many problems have arisen, particularly in connection with artifical 
satellites, for which this basic supposition is not valid and which require a more accurate treat- 
ment. The original approach is now of mainly historical interest: it permitted us to understand 
the basic features of the secular planetary perturbations of minor and major planets. 

The Lagrangian theory cannot fully solve the problem of the existence of the mean motion of 
the node and of the perihelion. Important theoretical progress in solving this problem has been 
achieved by Bohl (1909, 1912), Jessen (1935), Weyl (1938, 193 9) and Tornehave and Jessen (1945). 
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Skripnichenko (1968) applied their ideas to determine the mean motions of the ascending nodes of 
Venus, Earth, and Mars, and the mean motions of the perihelia of Venus, Earth, and Uranus,- 

Hagihara (1928) and Kozai (1954) both developed analytical theories of the secular perturba- 
tions of higher orders and higher degrees. They used Baker’s (1916) method of integrating the 
differential equations, and the classical expansion of the secular disturbing function in powers of 
the eccentricities and inclinations (which thus are assumed to be small). The second work by 
Kozai (1962) treats the secular perturbations of large eccentricities and inclinations under the 
assumption that the motion of the disturbed body is circular. 

With the advent of electronic computers these restrictions became unnecessary. By the use 
of step-by-step numerical integration, it became possible to penetrate the problem more deeply 
and to obtain a series of interesting disclosures about the secular behavior of the orbits of aster- 
oids and of satellites in cislunar space. Musen (1963) and Hamid* suggested the application of 
the Gauss (1818) method based on averaging the components of the disturbing force over the orbits 
of the disturbed and disturbing body. We have two basic modifications of Gauss’ method. The 
first was developed by Hill (1882) and re-discussed by Calladreau (1885). The second method 
was developed by Halphen (1888) and has been re-discussed by Goriachev (1937) and Musen (1963). 

Hamid* applied Hill’s method to compute the secular effects in the motion of comets. Hal- 
phen’ s method was used by Shute (1964) to compute the secular lunar effects in the motions of 
artificial satellites in elongated orbits in cislunar space, and by Smith (1964) for the motions of 
Enke’s comet and minor planets. Recently Musen and his associates have applied it to the inves- 
tigation of Hirayama families of minor planets. The results will be published in subsequent pa- 
pers. The author (1963) undertook the modernization of Halphen’s method, presenting it in terms 
of dyadics and vectors. In the present article we discuss Hill’s method in a similar manner. It 
is of interest to note that a careful re-examination of Hill’s original scaler development reveals 
its intimate connection with the vector algebra of Minkowski’s pseudo-euclidean three-dimen- 
sional space M 3 . To clarify the geometrical aspects of Hill’s theory we use in the present ex- 
position the algebra of dyadics and vectors in M 3 , and, when necessary, the algebra of dyadics 
and vectors in euclidean space E 3 . We propose a new symmetrical computational scheme in 
terms of traces of dyadics in M 3 . For the computation of elliptic integrals which appear in our 
exposition we propose Cody’s (1965) highly accurate approximations by Chebyshev polynomials. 

We have previously experienced difficulties in computing the secular effects for nearly cir- 
cular orbits when the eccentricity became negative and the Laplacian vector abruptly changed 
direction. A similar situation arose with small orbital inclinations. We develop in this article 
the differential equations for the secular perturbations in vectorial elements which are non- 
singular for small e and i. These elements are s = e//i - e 2 P 1 where P 1 is the unit vector di- 
rected toward the perihelion, in place of e and 77, and Gibbs’vector q (Gibbs, 1901) which defines 
the rotation of the orbital plane from its initial position to the position at a given time and which 
replaces the standard elements i, n, and a. The introduction of Gibbs’ vector (Musen, 1961) leads 


*Private Communication, 1963* 
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to differential equations for the perturbations of the position of the orbital plane in which the prob- 
lem of small divisors in form of the sine of the inclination is removed. The semi-major axis is 
not affected by the perturbations of zero rank; thus the computation of da/dt can serve as a check 
on the applicability of the method. 

We consider the case of only one disturbing planet; generalization to the case of a planetary 
system is not too difficult. The presented theory is applicable not only to the planetary case but 
also, for a more limited interval of time, to cometary orbits. It is applicable also to the orbits 
of space probes not approaching the Moon too closely. 


SOME BASIC RULES OF THE VECTOR ALGEBRA IN THE PSEUDO-EUCLIDEAN SPACE M 3 

For the sake of completeness of the exposition, we state briefly and without proof those rules 
of the vector algebra in m 3 which differ from the analogous rules in E 3 . We shall require these 
rules in our exposition of Hill T s theory. 

Let e x ,e 2 , e 3 be the basic unit vectors in M 3 . The fundamental multiplication table peculiar 
to M 3 (in fact, the conditions of pseudo- orthogonality) is 



e 2 • e„ = + 1, e 3 • c 3 = - 1, e. • e. = 0 (i 4 j) 

(1) 

? i x V + e 3 , 

? 2 x? 3 = -e i , e J x? i = - e 2 , e.xe.^0. 

(2) 

The dot product of two vectors 

® = a i *1 + a 2 ? 2 + a 3 ? 3' 



b = bj ?! + b 2 e 2 + b 3 e 3 , 


in m 3 , in agreement with (1), is given by 



a • b = a t bj + a 2 b 2 - a 3 b 3 ; 

(3) 


in particular 


— * — * 9 2 2 

a * a = + a^ - a‘ 


The cross product in M 3 , in accordance with (2), can be expressed as 


a x b = 


b i b 2 ~ b 3 
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The expansion formula for the double cross product in M 3 is 

ax (b x c) = (a ■ b) c - (a • c) b W) 

with the sign of the right side opposite to the sign in the analogous formula in e 3 . 

A linear vectorial transformation from the basic set ,e 2? e 3 ) to a new basic set (e^e^e*) 
which leaves (3) invariant is called a pseudo -orthogonal transformation. The idemfactor (the unit 
matrix) in M 3 has the form 


E = e. 


+ e 0 e, - e- e.. 


where the products of vectors are dyadic. From two representations of the idemfactor 


' 1 l + 2 " 3 c 3 “ 1 ^ 2 2 3 3 


and making use of (1) and of the analogous conditions for e*,e*, e*, we obtain the decompositions 


e . = a. j e t + a, a e 2 - a. 3 e 3 , 


e i _ a i i e l + a 2i e 2 a 3i e 3 1 


( 5 ) 

(6) 


where a. . = e* • e From (5) and (6) we have the formulas for the transformation of coordinates 

Col Of*, £*, -<f 3 ) = (0...) • col (£j, g 2 , £ 3 ), (7) 


col (£,. z 2 . -i 3 ) = (a..) • col (s'* , $ * , s" 3), 

(i, j = 1,2,3) 


(8) 


where i is the row index and j is the column index. The multiplication of matrices and column 
vectors in (7) and (8) is performed in the standard manner. The conditions of pseudo-orthogon- 
ality are 


x. , a., + a.„ a. « - a., a., = 5 . . , 
il 3 1 i2 j 2 i3 j3 1 j ’ 


. . a, . + cu . a_ . - a, . a, . = 5 . . 

1 1 1 j ~ 2 j 2 j 3 1 3 j 1 j 


( 9 ) 


where Kronecker deltas in M 3 are defined as 

+ 1 i ~ ) ~ 1 , 2 

h H = 0 1 ? 3 

- 1 i = j = 3 . 

Hill, in his exposition, resorts to (9). We do not use these relations in the present article. 
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The normal form of the dyadic 


<D = 2j bj + a 2 b 2 + 1 3 b 3 = 


7 c. . e. e. 

/ j 1 J 


(i, j = 1,2,3) 


in m 3 space is 


^ \ e l e i ~ S e 2 e 2 + ^3 e 3 e 3 ’ 

where k % , ^ 2 , \ 3 are the roots of the characteristic equation. In agreement with Callandreau we 
write this equation in the form 


< p (*0 = 


+ k 

C 1 2 

C 1 3 



21 

C 22 + * 

C 23 

= 0, 

(10) 

31 

C 32 

C 33 “ X 




which when expanded becomes 


<p(\) = \3 + *2 + f 2 ^ + f 3 = 0, 

where f x ? f 2 , f 3 are the scaler invariants of In M 3 they have the form 

fj = aj • bj + a 2 • b 2 + a 3 • b 3 = c n + c 22 - c 33 , 


(a, x a 2 ) • 

(b 1 x b 2 ) - (a 2 x a 3 ) 

• (b 2 x 

b 3 ) - (a 3 x a,) 

• (b 3 X bj) 

c u C 12 

C 22 C 2 3 

C 33 

C 3 1 


C 21 C 22 

C 32 C 33 

C 1 3 

C u 



C 11 C 1 2 C 1 3 

C 2 1 C 22 C 23 

C 31 C 32 C 33 


f 3 = - (a 1 ■ a 2 x a 3 ) (bj 


b 3 ) 


Cayley's identity 


O 3 - f x $ 2 + f 2 <J> - f 3 E = 0 
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remains valid in m 3 , as does the identity 

, E \ 2 - (<f -E f,) A. + ($ 2 - f, « +E f 2 ) (11) 

($ + X E )" 1 = , 

\ 3 + f x A 2 + f 2 K + f 3 

provided that A is not a root of the characteristic equation of <I>. 


COMPONENTS OF THE SECULAR DISTURBING FORCE 

To achieve compactness and symmetry of formulas, we use notations for the vectorial ele- 
ments of the disturbed and the disturbing body which differ from the standard notations. 
Substituting 

FT Aj (cos e' - e') + A' sin e' ( 12 ) 

into 

A = r' - r , 


we obtain 


A = A| cos e 1 + A' sin e* 




where p - r + e' is the position vector of the disturbed body relative to the center of the orbit 
of the disturbing body. We assume e' remains greater than zero. Decomposing p along the axes 
Pj' , P 2 , P 3 ', we obtain 

^ “ X 1 + X 2 ^2 + X 3 ^3 ’ 


where 


p ; 


+ e 


• P 2- 


P 3 ’ 


and (Hansen, 1957) 


A 2 =a' 2 (-y 0 + y 2 cos 2 e' - 2 y, cos e ' - 2 ft Q sin e ' ) , 


where 


,2 "i , 1 2 P 2 P x 2 71 - e ' 2 

>2 = e • :r — ■ >0 = 1 ~ e + — ■ A = : 


In the actual computations and programming, the system of formulas 


y 2 = e' 2 , y x - e # + a — • Pj, y Q =- 1 - 2e' 2 + a 2 — + 2e' y y ^-al - P 2 /l-e' 2 


(13) 


(14) 
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where 


is preferable to (13). The expressions 

1= 1 - e cos e, — = P, (cos e - e) + P 0 /l - e 2 sin e 
a a 1 2 

are to be substituted into (14). 

In the ensuing exposition it will be convenient to express A 2 in terms of the position vector 

f = + f 2 ? 2+ ^3 ? 3 ^ 


of the point 


= r cos e' , = r sin e', g 3 = r 


(16) 


lying on the surface 


£•£=£?+ f J - f ? - 0 


in the space M 3 ; here r is a factor to be discussed later. 
Introducing the M 3 dyadic 


<t = y 2 ?i + y, (e, e 3 + e 3 ?,) + /?„ (e 2 e 3 + e 3 e 2 ) + y 0 e 3 e 3 


and taking (1), (15), and (16) into account, we obtain 




(17) 


(19) 


By means of a pseudo- orthogonal transformation we can reduce (17) to the normal form 

<D = G, ej e* - G 2 e * + G 3 e;e^, 

where -Gj , +G 2 and +g 3 are the roots of the characteristic equation 


<p ( K ) 


+ y 2 0 - y\ 

o k -P 0 

yi Po x - 


= 0. 


(19) 


( 20 ) 
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To help the reader compare the original expositions of Hill and Callandreau with the present ex- 
position, we are using Hill T s notations for the characteristic roots. In the expanded form we have 

<p 00 = [*■(*•- Vo) ( k + y 2 ) + r? ^ = o, (21) 

or 

cp (\) = + fj \ 2 + f 2 A. + f 3 = 0 , (22) 


where 

f l=r 2 -r 0 > f 2 = ^1 +^0 - r 0 7 2 ’ f 3 =/ 3 0>2 

are the scalar invariants of <P. 

Hill has shown that all roots of (22) are real, 
cations. From (21) we obtain 

<p (-r 2 ) = -r? y 2 . <p (0) = + /3q r 2 . 

and, taking (13) into account, 



a' 2 


We repeat his proof here, with slight modifi- 

( + r 0 ) = ^0 fro +> , 2 > + ^>V 


Thus 


<P (-7 2 ) < 0, cp (0) >0, cp (1 - e' 2 ) <0, cp (+7 0 ) > 0. 

Consequently, the characteristic equation has three real roots, - Gj , +G 2 , +g 3 , located in the intervals 

(-7 2 . 0), (0, 1 - e' 2 ) and (1 - e' 2 , y 0 ) • 

The root -G, can be easily obtained by the method of iteration. From (22) we have 

Q, = li 

<? - f, O, * f, 

and, because of the smallness of y 2 , we obtain G : in the case of a minor planet, or a satellite in 
cislunar space, after only a few iterations. After G 1 is computed, the remaining roots + g 2 and 
+ g 3 can be determined from the reduced equation 

\ 2 + (f t - Gj) A. + (G 2 - fj Gj + f 2 ) = 0. 

In fact, (22) appears in Hansen's work (1857) in connection with the numerical expansion of the 
disturbing function, where the use of the iteration procedure is also suggested. 
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In our exposition we shall need only the roots of (22), so we can dispense with the determina- 
tion of the matrix of the pseudo -orthogonal transformation, as well as with the actual determination 
of the factor t. We select r such that the new coordinate g* of the point (16) will be equal to unity. 
Then, taking £ = 0 into account, we can set 

= cos T, f * = sin T, f*=+l. ( 23 ) 

By substituting 

£ = e * cos T + e* sin T + e* ( 24 ) 

and (19) into (18) we obtain 

r 2 ^ = Gj cos 2 T - G 2 sin 2 T + G 3 . 

Differentiating the identity 

e* cosT + e * sinT + e 3 = r (e t cos a 1 + e 2 sin a 1 + e 3 ) , ( 24 f ) 

and taking (2) into account, we deduce 

£ x e* dT = - dr + £x e 3 d e ' 


or 


?X <?x SJ) dT^x (? X e 3 ) d€ # , 


and, after the application of (4), 


d T = rde'. 

The disturbing force averaged over the orbit of the disturbing body is 

*2rr 

. i i 

F n = 


. 


Taking (12) into account we obtain 


where 


r 2 

?.a» = _L f 
0 ^ 1 


(Nj cos a 1 + N 2 sin a* + N 3 ) 

d g' . 

A 3 


Nj = Pj • u°, N 2 = P' • u° /l - e' 2 , N 3 = - Pj • u° e' - a i_ • u°, 


(25) 


(26) 
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for the component of F 0 in the direction of a unit vector u° . Then, taking (18) and 

£ 1 - e # sin e* = g', — = 1 - e' cos c' 

a' 

into account we have, from (26) 


where 


a ' 2 F 0 




a • * • ?) 3/2 


de\ 


n = - ( N J + N 2 ? 2 - N 3 ?3> (e 3 + e' ej) 


( 27 ) 


(28) 


is an M 3 -dyadic. The trace of n is 


Pl = - (N 3 + e f N x ). 


From (7), (8) and (23) we obtain 


- 1 = r (o- 31 cos e* + a 32 sin e* + a 33 ), 

- r = a u cos T + a 23 sin T + a 33 . 


To avoid a contradiction between these last two relations we must conclude that r cannot become 
zero or infinity, but oscillates between two fixed limits which, with the proper choice of the "’di- 
rection cosines”, can be assumed to be positive. From (25) it is evident that T is a monotonically 
increasing function of e ' and that £ is a periodic function of T with period 2 tt-. In addition, it is 
clear from (24 r ) that when t ' covers a full period T also covers a full period. 


When (25) is taken into account, (27) becomes 


J r 2 
o 


f • n * £ 


■ dT, 


(29) 


(£ * <J> * O 372 

where $ is taken in its normal form (19), and correspondingly n is taken in the form: 


n 


— ► * — * * i—i — ** — * * i—i — ► * — * + \ ■ — * * — * * 

1 e l £ 1 + r 2 e 2 e 2 + r 3 e 3 e 3 + A ii e i e j’ 

with ~i in the form (24). But, since 


(30) 


i, j = 1 

i =* J 



sin T cos T dT 

(G t cos 2 T - G 2 sin 2 T + G 3 ) 3/2 


cos T dT 

(Gj cos 2 T - G 2 sin 2 T + G 3 ) 3/2 


= 0, 


= 0, 
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and 


1 r 27r 

2 71 Jo (Gj CO 


sin T dT 


= 0, 


s 2 T - G, sin 2 T + G,) 3/2 


all terms having coefficients A ;j disappear from (29) and we have simply 


a ' 2 f • n° = J_ f 

0 2 - Jo (Gj 


27T cos 2 T + r 2 sin 2 t + r 3 


T - G 2 sin 2 T + G-) 


3/2 


dT, 


or, taking into account the symmetry of the integrand, 

r n/2 r x cos 2 T + r 2 sin 2 T + P 3 


■' 2 v=°4 J 


Jo (Gj cos 2 T - G 2 sin 2 T + G 3 ) 3/2 

Consequently, we have to retain only the purely quadratic portion 


.dT. 


r — * * — * * i-. — * * — * * i — i 

1 e l e i + r 2 e 2 e 2 + r 3 


(31) 


in the transformed n. From 


<P + \E = (G 1 + k) e J e* 


“ (G 2 - k) e* e * + (G 3 -\)P;e* 


we obtain 


(<D + \ E) 


.1 = g l g l g 2 ^2 g 3 g 3 

G 1 + A. G 2 - A. G 3 - A. 


provided that A. is not a root of (22). From this last equation and from (30) we have 

r. 


n • (<D + A. E) 


p r 

2 1 3 

e „ e n - — e , e _ + • • * 


J 1 _r " 2 

and, forming the trace of the last expression, 


_ i _ i —**—*+ 

~ gTTX e ’ e * "gT^I " 2 " 2 "gT^I " 3 " 3 


{n • (<j> + xe )- 1 } 1 = 


r i 

Gj + \ 




(32) 


We next obtain this trace using the unreduced forms of <1> (17) and n (28) and the identity (11); 
and then, by comparing the two expressions for the trace, we determine r x , r r From (17) we 
have 


= (y\ - 7?) e, e„ - y t /3 0 (e t e 2 + e 2 + y 1 (y 2 - y 0 ) (e 1 e 3 + e 3 e t ) 
” Po e 2 e 2 “ Po 7q ( e 2 e 3 + e 3 e 2^ + (P 0 * 7 i ” 7 q) e 3 e 3 » 


and, from this last equation, and (17) and (28), we obtain 

n • <t> = (N t + N 2 e 2 - N 3 ? 3 ) (S u e 3 + s 12 e 2 + s 13 e 3 ) , 


(33) 
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where 


s n = y* - e '?v s i2 =^0 > s i 3 = >0 - e 'v 

and 

n • «I> 2 = (Nj + N 2 e 2 - N 3 ? 3 ) (s 21 ?j + s 22 e 2 + s 23 e 3 ), (34) 

where 

s 2i = y\ (y 2 -y<> > - e ' (?! -y?). s 2 2 = -£ 0 + e ' ?iA>> s 2 3 = (&l + y\ ~y\ ) - e ' n (r 2 -? 0 )• 

The traces s 3 and s 2 of (33) and (34) are 

S 1 = (n • ®)l = N 1 s n + N 2 S 12 + N 3 S 13 - 

s 2 = <n • « >2 ) ] = Nj s 23 + n 2 s 22 + n 3 s 23 . 


Multiplying n by (11) and taking the traces of the left and right sides we deduce 

<n - 11 - fj) 

A. 3 + f t \ 2 + f 2 X + f 3 

Then by equating the right sides of (32) and (35) and making use of l’Hopital’s rule we obtain: 


(35) 


p, G 3 + ( S j - f 3 p,) G t + (s 2 - s, fj + p t f 2 ) 
3 G 2 - 2 f, Gj + f 2 

Pi G 2 ~ < S 1 - f l Pi) G 2 + ( S 2 ~ S 1 f l + Pi f 2 ) 
3 G 2 + 2 fj G 2 + f 2 


r 3 = 


p, G 2 - (s 3 - f 3 Pl ) G 3 + (s 2 - Sj fj + Pl f 2 ) 
3 G 2 + 2 f j G 3 + f 2 


With fj , r 2 , r 3 now determined, we can rewrite (31) in the form 


a' 2 F Q ■ ?° (G, + G 3 ) 3/2 =1 f 

Jo 


2 C n/2 r t cos 2 T + r 2 sin 2 T + ^ 


(1 - k 2 sin 2 T) 3/2 


dT, 


(36) 
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where 


k 2 = 


G 1 + °2 
G 1 + G 3 


Making use of the formulas 



dT _ E (k) 

(1 -k 2 sin 2 T) 3/2 k' 2 



sin 2 T dT 
(1 - k 2 sin 2 T) 3/2 

cos 2 T dT 
(1 - k 2 sin 2 T) 3/2 


E (k) - k' 2 K (k) 
k 2 k' 2 


_ K(k) - E(k) 
k 2 


where k' 2 = i - k 2 and K(k) and E(k) are the normal elliptic integrals of the first and second kind 
respectively, we can put (36) into a form convenient for the numerical computations: 


7T a' 2 (G t + G 3 ) 


3/2 


( r 2 + r 3 )-^ + 

k 


f 2 ) 


K - E 
k 2 


For the numerical evaluation of elliptic integrals, we recommend their highly accurate approxi- 
mations of Chebyshev polynomials as obtained by Cody (1965). 

For the special cases of (26), 


J f%2 

0 


cos e - e 


A 3 


dg' 


K 2 = 


1 f 2n sin e 1 

277 Jo A 3 


dg' 


J f2n 

o A 3 


we obtain 


K i = ■ 


V a' 2 (G 1 +G 3 ) 3 


( r i 5 
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II II III 


where 


(s 21 e s 23 ) + (s lx e s 13 ) (G x - f x ) 

r n = + 

3 Gj - 2 fj G 1+ f 2 

(s 2 i - e s 23 ) - (s n - e s 13 ) (G 2 + fj) 

r i 2 = + 

3 G 2 2 + 2 fj G 2 + f 2 

r = - (s 2i ~ e ' s 23 ) - ( s n - e ' s i 3 ) (° 3 + f l> . 

3 01 + 2 ^ 0 ,+^ 

and 


a ' 2 (G, + G 3 ) 3 


^22 + ^ 23 ^ — — + (^21 

k * 


K - E 
h 2 


where 


S 22 + s i 2 ( G i ^ 1^ 

3G?-2f 1 G 1 + f~ ’ 


S 22 ~ 5 1 2 (°2 + f l> 

3 G 2 2 + 2 f t G 2 + f 2 ’ 


S 22 S 12 ^3 + ^l) 

3 G f + 2 fl g 3 + f a ’ 


and 


K 


3 


2 


(G t 


+ G 3) 3 



^ 33 ) — + (r 31 

33 k , 2 


^ 32 ^ 


K - E 
k 2 


where 


31 


G 1 “ ( S 13 + ^1) G 1 ” ( S 23 " S 1 3 “ ^2) 

3 G 2 - 2 fj G x + f 2 
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II 


p _ ^2 + ( S 13 + ^2 “ ( S 23 “ S 13 ^1 “ ^ 2 ) 

3G 2+2f 1 G 2 +f 2 

p _ + ^3 + ( S 13 + ^l) ^3 “ ( S 23 " S 13 ^1 “ ^ 2 ) 

3 G 3 2 + 2 f x G 3 + f 2 

EQUATIONS FOR SECULAR VARIATION OF ELEMENTS 

To avoid the difficulty associated with division by a small eccentricity, when the elements n 
and e are used separately, we suggest instead the use of the vectorial element 


s 



which is intimately related to Hamilton's vector and to the element (h/h 0 ) ePj of Hansen's lunar 
theory. The projections of s on the axes of the ideal reference frame (a frame rigidly connected 
with the osculating orbit plane of the disturbed body) are 


s 



cos y, 



sin 




where y is the true orbital longitude of perihelion of the disturbed body reckoned from the depar- 
ture point of the ideal system. 

The differential equation of motion of the disturbed body, as referred to the ideal system, can 
be written in the form 


d 2 r 


dt 2 


}jl ~r _ fj. m' 
r 3 M + m 


(F), 


where (f) is the projection of the total disturbing acceleration 


( 37 ) 


F = 


A 3 r ' 3 
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onto the osculating orbit plane. For our purpose it is convenient to write the Laplace qua si -integral 
of (37) in the form 


R„ 



= 0 . 


(38) 


The unit vector R 3> normal to the orbital plane, is a constant in the ideal system. Differentiating 
(38), taking into account (37) and 


d 1 m / /yLi^'*^3 xr d-r° n° v'fj. p 

d t v/ - M + m p dt r 2 

and considering the fact that a is not affected by the perturbations of zero-rank, we average the 
result over the orbits of both bodies and obtain for the perturbations of zero-rank in s : 


d s _ m' n a 2 1 
d t M + m 2 7T 



K • F.dg. 


(39) 


where 


5 i 5 2 - g 2 5, + 


and 


f o = a ; k i 


+ A' K 2 - r K 3 . 


(40) 


In the frame of the present theory the ideal system (r 1? r 2 , r 3 ) is affected only by the pertur- 
bations of zero-rank. The averaging process in (3 9) can be greatly facilitated by replacing the 
integration with respect to g by integration with respect to the true orbital longitude A. Making use 
of the geometric relations 


*1 = — /TT T 2 . !* = o, 

r 2 9 g 


where v is the true anomaly, 


r = r (R t cos A + R 2 sin A) , 


- R, 


sin A. + R 2 cos K, 


(41) 
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( 42 ) 


we obtain, instead of (39), 


d s _ m # n a 2 1 
d t M + m 2 77 


f 

Jo 


(Wj + W 2 K 2 + W 3 K 3 ) d\, 


where we set 


W. = ( I ) 2 V • AJ /l + S 2 , (j = 1, 2) 


w 3 = - ^ ^ v • 7 /i + s 2 . 


In scalar form, noting that in the ideal system R i (i = 1, 2) are constant, we have 


d s i 

dt 


m' n a 2 1 
M + n 2 77 



Sij K j d ^’ 


(i = 1. 2) 


(43) 


where 


A' * R 9 + — (1 + s 2 )n°cos \ ( L\ /l 

J L a J \ a / 


+ s" 


(j = 1. 2) 


S 2j - A! * -Rj + ^ (1 + s 2 )n°sin X.J /l + s 2 , 


'13 




/l + s 2 s in \ , 


(j = 1. 2) 


S 23 = + a /l + s 2 cos 
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We determine the perturbations of zero-rank in the position of the orbital plane by means of 
Gibbs’ vector 

3 

i = 1 

decomposed along the s. (i = 1, 2, 3); s i are the initial values of the unit vectors of the ideal refer- 
ence frame. The matrix of rotation A of the ideal frame from its initial position to its position at 

the time t, as expressed in terms of q, has the form (Gibbs, 1901) 

A = I + — [q x I + q x (q x I) ] = I + — (q x I + q q) . 

1 + q 2 1 + q 2 1 + q 2 

In addition, we shall use the dyadic 

H = i(l +q 2 ) (I +A) = I + qx I + qq, 

in the matrix form 

== 2 Z (ii j = i - 2 ’ 3) - 

i , j 

where 

^ + q ? ’ ^ 1 2 = ^ 1 ^2 — ^3 ’ ^13 = ^ 1 ^3 ^2 ’ 

^21 “ ^1 ^2 + q 3 ’ ^22 = ^ + ^2 * ^23 “ ^2 ^3 ” q 1 ’ 

^31 =%% -<* 2 ' «32 = ^3 + % • ^33 = 1+C lI- 

The unit vectors Rj (i = 1, 2, 3) are given by the formula 

R. = S. + — - — [q X S. + q X (q X S.)]= 1 ~ q S. + — (q x S. + q q • S ; ) . 

1 + q 2 1 + q 2 1 + q 2 
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If we set 


A ; i=A ; • r,, 


it now follows from (44) that 


A ii = 4 X I ■ + ^ x a; + Qi g • a.). 

1 + 1 + q z 


(45) 


After these A/, are determined using (45), we compute s L . (i, j = 1, 2) by means of the formulas 


S Xj = A! 2 + — (1 + s 2 ) (- A! x sin k + Aj 2 cos k) cos 


S 2j = -A' x + — (1 + s 2 ) (- A' l sin k + A' 2 cos A.) sin A.J 



/l + s 2 , 



These expressions are to be substituted into equations (43). 

We have previously established (Musen, 1961) the following differential equation for the per- 
turbations in Gibbs’ vector: 


dq _ 1 

TT" 2 


M + m /j e 2 a 


r ( F • r 3 ) = 


(Sj cos A. + S 2 sin A.) . 


Taking the average over the orbits of both bodies and replacing the integration with respect to g 
by integration with respect to the true orbital longitude, we obtain, for the zero-rank perturbations 
in q 


dq 

dt 


1 m' n a 2 

2 M + m 



( a ; 3 k i 


+ A 23 K 2 ) 5 ■ (Sj cos k + S 2 sin A) d A. 


(46) 

By projecting (46) on Sj (i = 1, 2, 3) we obtain the scalar equations to be used in the actual compu- 
tations: 


dc h 

dt 


M + m 


•(1 + s 2 ) 


Jlrr 

1 f /r_ \ 3 

2 7T l a ) 

Jo ' / 


( a ; 3 r i 


a ' 3 k 2 ) (g n cos \ + sin A) d 


A.. 
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The components of s in the inertial system of coordinates and, if necessary, the standard elliptic 
elements e, i, 0, can be easily determined from the equation 


S = A • (Sj + s 2 s 2 ) . 


As a check we can compute the derivative of the semi-major axis. Its smallness guarantees 
the accuracy of the theory. From 


d a 

cTT 


2 m' n a 3 1 
M + m 2 7T 



(r 0 e sin e + n 0 /l - e^) 


de, 


and taking into account (40) and (41) and 


d 


r d k 
a /TT72 ' 


we obtain 


d a 
d t 


(“l K 1 + a 2 ^2 + a 3 K 3>^, 

0 


where 


a . 

) 


: i(/r 


+ s e sin g 


r° • a; + f?° - a:), 


(j = 1. 2) 


a 


3 


a 



e sin e , 


and 


a; = + a' 


cos \ + A'. 2 sin A , 


A. 


A^ sin K + A' 2 cos A, 
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r/a is given in terms of s by the equation (42) and we have also 


r / 

e sin e = — (s. 

a 1 


sin k - s 2 cos K) . 


CONCLUSION 

We have developed a set of symmetric formulas for computing the zero-rank perturbations 
in the framework of Hill's theory. This symmetry facilitates the optimization of programming. 
The vectorial elements of motion are chosen such that the system can also be used when the 
eccentricity or the inclination of the orbit becomes small but oscillates in a wide interval. 

Hill's and Halphen’s methods both become inapplicable if two orbits come too close together. 
In Hill’s method the numerical difficulty caused by the proximity of orbits appears as a small 
numerical divisor k'. Thus, in Hill’s method, the difficulty can be watched more easily and di- 
rectly than in Halphen’s method. Yet we hesitate to give a definite preference to either, partly 
because we have succeeded in applying Halphen’s method to determine the long range effects in 
the orbits of planets, comets, and space probes. However, Hill’s method is appealing to the ce- 
lestial mechanician because of its geometrical simplicity and elegance. Almost every transfor- 
mation in Hill’s method has a direct geometrical or kinematical meaning. 

The methods of numerical averaging presently have certain advantages over purely analytical 
methods. They can treat a large range of eccentricities and orbital inclinations. They can also 
treat the free secular oscillations as well as the forced ones, and also their mutual cross- 
effects. At the present time, no analytical theory can do this completely. Perhaps this is because 
we are not in possession, and hopefully never will be (for the sake of science), of a general analy- 
tical solution of the many body problem. 
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